1 Goal / questions

  1. What probability distribution should we use?
  2. How should we estimate parameters? Is the variation different between rivers?
  3. Do we need to take temporal correlation into account?

2 Counters

2.1 Data

library(tidyverse)
library(modelr)
library(zoo)
library(patchwork)
library(mar)

mar <- connect_mar()
default_theme <- theme_set(theme_bw())

gogn <- read_delim("teljaragogn.csv") |>
  filter(!is.na(fjoldi)) |>
  filter(!vatnsfall %in% c("Selá","Laugardalsá")) |>
  #filter(!(ar < 1992 & vatnsfall=="Blanda")) |>
  mutate(ln = log(fjoldi)) |>
  group_by(vatnsfall) |>
  mutate( m5 = lag(rollmean(ln, 5,fill=NA, align = "right"))) |>
  ungroup()
gogn |>
  ggplot(aes(ar, fjoldi, col=vatnsfall)) +
  geom_point() +
  geom_path() +
  facet_wrap(~vatnsfall, scales = "free") +
  theme(axis.text.x=element_text(angle=45, hjust=1)) 

gogn |>
  ggplot(aes(ar, fjoldi, col=vatnsfall)) +
  geom_point() +
  geom_path() +
  facet_wrap(~vatnsfall, scales = "free") +
  theme(axis.text.x=element_text(angle=45, hjust=1)) +
  scale_y_log10()

gogn |>
  group_by(vatnsfall) |>
  summarise(start=min(ar), stop=max(ar), fjoldi_ara=n(), medal_fjoldi=mean(fjoldi), log_m_trans=exp(mean(ln)), log_medal= mean(ln), stadalfravik= sd(fjoldi), log_stadal=sd(ln)) |>
  arrange(medal_fjoldi) |>
  kableExtra::kbl(digits = 2) |>
  kableExtra::kable_styling("striped", full_width = F)
vatnsfall start stop fjoldi_ara medal_fjoldi log_m_trans log_medal stadalfravik log_stadal
Krossá 1998 2023 25 226.92 186.54 5.23 141.18 0.66
Vesturdalsá 1994 2023 25 229.72 208.91 5.34 95.63 0.47
Búðardalsá 2016 2024 9 292.33 251.35 5.53 163.67 0.60
LangáSveðja 2000 2023 24 962.71 861.83 6.76 441.17 0.49
Blanda 1982 2023 42 1518.60 1174.30 7.07 1059.22 0.78
NorðuráGlanni 2002 2024 23 2188.57 1969.47 7.59 1035.97 0.48
Elliðaár 1933 2023 91 2378.84 2088.49 7.64 1202.84 0.53 
gogn |>
  group_by(vatnsfall) |>
  summarise(start=min(ar), stop=max(ar), fjoldi_ara=n(), medal_fjoldi=mean(fjoldi), log_m_trans=exp(mean(ln)), log_medal= mean(ln), stadalfravik= sd(fjoldi), log_stadal=sd(ln)) |>
  ggplot(aes(log_medal, log_stadal))+
  geom_point() +
  geom_smooth(method = "lm")

3 Fish counters - distribution

mNorm <- lm(fjoldi~vatnsfall, data=gogn) #Meðaltöl í fjölda
mLNorm <- lm(ln~vatnsfall, data=gogn) #Meðaltöl á log
mLag <- lm(ln~vatnsfall + m5, data=gogn)
mLag2 <- lm(ln~vatnsfall + vatnsfall*m5, data=gogn)

gogn |>
  spread_residuals(mNorm,mLNorm, mLag, mLag2) ->
  gogn

3.1 Normal distribution

plot(mNorm)

shapiro.test(gogn$mNorm)
## 
##  Shapiro-Wilk normality test
## 
## data:  gogn$mNorm
## W = 0.91091, p-value = 9.62e-11

3.2 Log-normal

plot(mLNorm)

shapiro.test(gogn$mLNorm)
## 
##  Shapiro-Wilk normality test
## 
## data:  gogn$mLNorm
## W = 0.98565, p-value = 0.01663

4 Fisktalning - Minni

mLNorm <- lm(ln~vatnsfall, data=filter(gogn,!is.na(m5)))
anova(mLNorm, mLag)
anova(mLag, mLag2)

4.1 Partial autocorrelation of residuals from ANOVA

for (v in unique(gogn$vatnsfall)){
  #acf(gogn$mLNorm[gogn$vatnsfall==v], main=v)
  pacf(gogn$mLNorm[gogn$vatnsfall==v], main=v)
}

4.2 Partial autocorrelation of residuals from ANCOVA (5-year running mean)

for (v in unique(gogn$vatnsfall)){
  #acf(na.omit(gogn$mLag[gogn$vatnsfall==v]), main=v)
  pacf(na.omit(gogn$mLag[gogn$vatnsfall==v]), main=v)
}

4.3 Standard deviation

summary(mLNorm)$sigma
## [1] 0.5679721
summary(mLag)$sigma
## [1] 0.5038929

5 Salmon catches

5.1 Data

listi <- c("Andakílsá"=12,"Blanda"=158,"Elliðaár"=411,"Norðurá"=1861, "Krossá"=1462,
           "Langá"=1568,"Miðfjarðará"=1733, "Hítará"=956,"Sog"=2271,"Stóra-Laxá"=2339,
           "Ásgarðslækur"=284,"Skaftá"=2164,"Héraðsvötn"=955,
           "LaxáLeir" = 1589,"LaxáAðal"=1591,"Þverá"=2870,"LaxáKjós"=1590, "Álftá"=69,
           "Fnjóská"=537,"Grímsá"=749,"Haffjarðará"=835,
           "Hafralónsá"=851,"Haukadalsa"=883,"Hofsa"=996,"Laugardalsa"=1579,
           "Svalbarðsa"=2424, "Brynjudalsá"=257,"FáskrúðDölum"=465,"Hrútafjarðará" =1106,
           "HölknáÞist"=1216, "LaxáDölum"=1600,"Leirvogsá"=1615, "SeláVopn"=2080,
           "StaðaráSt"=2282)
veidibok_veidibok(mar) |> 
  left_join(tbl_mar(mar,"ops$johannes.VEIDIBOK_TO_WATER"), by=c("vatnsfall_id"="id")) |>
  left_join({fv_vatn(mar) |>
      select(vatn_nr, adalvatnsfall_nr, vatnsfall_til_nr)}, by=c("water_no"="vatn_nr")) |>
  left_join({fv_vatn(mar) |> 
      select(vatn_nafn, vatn_nr) |> 
      rename(adalvatnsfall = vatn_nafn)}, by= c("adalvatnsfall_nr"="vatn_nr")) |>
  filter(fisktegund_id == 91, ar <2023, ar > 1973) |>
  filter(adalvatnsfall_nr %in% listi |
           water_no %in% listi) |>
  filter(!(water_no==1568 & ar %in% 1989:1990),
         !(adalvatnsfall_nr==996 & ar==2008),
         vatnsfall!="Blöndulón",
         !(water_no==2574 & ar %in% 2008:2009), #Tungulækur hafbeit
         water_no!=922) |> #Hellisá flutningur á fiski
  #replace_na(list(sleppt_kodi=0)) |>
  group_by(adalvatnsfall_nr) |>
  collect(n=Inf) |>
  mutate(vatnsfall= str_flatten(sort(unique(vatnsfall)), collapse = " ")) |> 
  ungroup() |>
  count(vatnsfall, adalvatnsfall_nr, ar) |> 
  mutate(ln=log(n)) |>
  arrange(adalvatnsfall_nr, ar) |>
  group_by(adalvatnsfall_nr) |>
  collect(n=Inf) |>
  mutate( mInf = mean(ln,na.rm = T),
          m5 = lag(rollmean(ln, 5,fill=NA, align = "right")),
          m10 = lag(rollmean(ln, 10,fill=NA, align = "right"))) |>
  ungroup() ->
  test
test |>
  ggplot(aes(ar,n, col=vatnsfall)) +
  geom_path() +
  #geom_point() +
  #geom_path(aes(y=m5), linetype = 2) +
  #geom_path(aes(y=m10), linetype = 3) +
  facet_wrap(~adalvatnsfall_nr, scales="free_y") +
  theme(legend.position = "none")

test |>
  ggplot(aes(ar,n, col=vatnsfall)) +
  geom_path() +
  #geom_point() +
  #geom_path(aes(y=m5), linetype = 2) +
  #geom_path(aes(y=m10), linetype = 3) +
  facet_wrap(~adalvatnsfall_nr, scales="free_y") +
  theme(legend.position = "none") +
  scale_y_log10()

test |>
  group_by(adalvatnsfall_nr,vatnsfall) |>
  summarise(start=min(ar), stop=max(ar), fjoldi_ara=n(), medal_fjoldi=mean(n), log_m_trans=exp(mean(ln)), log_medal= mean(ln), stadalfravik= sd(n), log_stadal=sd(ln)) |>
  #arrange(medal_fjoldi) |>
  kableExtra::kbl(digits = 2) |>
  kableExtra::kable_styling("striped", full_width = F)
adalvatnsfall_nr vatnsfall start stop fjoldi_ara medal_fjoldi log_m_trans log_medal stadalfravik log_stadal
12 Andakílsá 1975 2022 43 213.98 168.78 5.13 177.14 0.65
69 Álftá á Mýrum Veitá 1974 2022 49 286.39 264.62 5.58 120.30 0.40
158 Blanda Svartá 1974 2022 49 1423.14 1197.31 7.09 923.57 0.60
257 Brynjudalsá 1975 2022 45 149.24 114.78 4.74 108.85 0.81
411 Elliðaár Elliðavatn 1974 2022 49 1055.63 967.19 6.87 430.47 0.44
465 Fáskrúð 1974 2022 46 225.67 206.44 5.33 100.18 0.43
537 Bakkaá Fnjóskadal Fnjóská 1974 2022 48 285.75 234.56 5.46 191.25 0.64
749 Grímsá Borgarfirði Tunguá Lundarreykjadal 1974 2022 49 1261.00 1178.09 7.07 445.98 0.39
835 Haffjarðará Hlíðarvatn Mýrum 1974 2022 49 926.22 849.43 6.74 423.62 0.41
851 Hafralónsá Kverká 1974 2022 47 280.28 244.80 5.50 120.35 0.62
883 Haukadalsá Haukadalsá Dýrafirði Haukadalsvatn 1974 2022 49 676.39 632.14 6.45 247.04 0.38
955 Héraðsvötn Hofsá Vesturdal Húseyjarkvísl Miklavatn Norðurá Skagafirði Sæmundará 1974 2022 49 260.96 208.03 5.34 170.40 0.71
956 Hítará 1974 2022 49 479.78 417.89 6.04 266.10 0.53
996 Hofsá í Vopnafirði Sunnudalsá 1974 2022 48 1004.00 844.16 6.74 539.68 0.65
1106 Hrútafjarðará og Síká 1976 2022 46 343.22 310.79 5.74 160.61 0.45
1216 Hölkná 1974 2022 48 95.38 81.28 4.40 47.58 0.63
1462 Krossá 1974 2022 49 117.43 96.62 4.57 73.69 0.65
1568 Langá 1974 2022 47 1540.55 1425.55 7.26 610.30 0.40
1579 Laugardalsá Laugardalsvatn 1974 2022 48 296.10 255.06 5.54 162.85 0.56
1589 Eyrarvatn, Þórisst. Geitabergs Laxá í Leirársveit Selós Svínadal Þórisstaðavatn Þverá Svínadal 1974 2022 47 965.79 907.32 6.81 341.87 0.36
1590 Bugða Laxá í Kjós Meðalfellsvatn 1974 2022 49 1334.08 1238.45 7.12 540.76 0.40
1591 Kráká Laxá í Aðaldal Mýrarkvísl Reykjadalsá og Eyvindarlækur 1974 2022 49 1786.92 1547.79 7.34 998.05 0.55
1600 Laxá í Dölum 1974 2022 47 1014.00 898.18 6.80 481.41 0.52
1615 Leirvogsá 1974 2022 49 461.37 414.10 6.03 214.18 0.49
1733 Austurá Miðfjarðará Núpsá Vesturá 1974 2022 49 1730.14 1453.36 7.28 1131.13 0.59
1861 Norðurá Borg 1974 2022 49 1667.98 1550.21 7.35 643.39 0.40
2080 Selá í Vopnafirði 1974 2022 48 1152.52 979.02 6.89 574.35 0.66
2164 Geirlandsá Grenlækur Holtsá Hörgsá Hörgárdal Hörgsá Síðu Jónskvísl og Sýrlækur Skaftá Tungulækur 1974 2022 49 62.41 55.75 4.02 29.91 0.51
2282 Staðará Steingrímsf 1974 2022 42 79.52 67.99 4.22 44.12 0.59
2424 Svalbarðsá Svalbarðsá Þistilsfirði 1974 2022 48 254.40 208.90 5.34 148.33 0.70
2870 Kjarrá Litla-Þverá Þverá Þverá Borgarfirði eystri 1974 2022 48 1927.25 1805.30 7.50 735.15 0.36
2959 Ásgarðslækur Búrfellslækur Sog Stóra-Laxá 1974 2022 48 844.92 775.01 6.65 405.08 0.41 
test |>
  group_by(vatnsfall) |>
  summarise(start=min(ar), stop=max(ar), fjoldi_ara=n(), medal_fjoldi=mean(n), log_m_trans=exp(mean(ln)), log_medal= mean(ln), stadalfravik= sd(n), log_stadal=sd(ln)) |>
  ggplot(aes(log_medal, log_stadal))+
  geom_point() +
  geom_smooth(method = "lm",se=FALSE) +
  geom_smooth(col="black", se=FALSE)

6 Catch data distribution

Skip checking the normal distribution (it’s obviously does not fit)

6.1 log-normal

#Líkön
mL <- lm(ln~vatnsfall, (test))
mLm5 <- lm(ln~vatnsfall + m5, data = (test))
mLm10 <- lm(ln~vatnsfall + m10, data = test)
mLNA <- lm(ln~vatnsfall, na.omit(test))
mLm5NA <- lm(ln~vatnsfall + m5, data = na.omit(test))
mLm5m10 <- lm(ln~vatnsfall + m5+m10, data = test)
mL0m5 <- lm(ln~m5, data = na.omit(test))
mL0m5m10 <- lm(ln~m5 + m10, data = na.omit(test))
mL0m10 <- lm(ln~m10, data = na.omit(test))
mLm5NAspec <- lm(ln~vatnsfall + vatnsfall*m5, data = na.omit(test))

test |>
  spread_residuals(mL, mLm5,mLm10,mLm5m10) ->
  test
plot(mL)

shapiro.test(test$mL)
## 
##  Shapiro-Wilk normality test
## 
## data:  test$mL
## W = 0.98438, p-value = 8.127e-12
shapiro.test(test$mLm5m10)
## 
##  Shapiro-Wilk normality test
## 
## data:  test$mLm5m10
## W = 0.99456, p-value = 0.0002248

Feitir halar, veiðihlutfall og léleg skráning?

7 Catch data lag 5 year, 10 year

drop1(mLm5m10, test="F")
drop1(mLm5NAspec, test = "F")

5 year mean with highest F-value

7.1 Standard deviation

summary(mLm5m10)$sigma
## [1] 0.4694884
summary(mLm5)$sigma
## [1] 0.4922756
summary(mL)$sigma
## [1] 0.5391005

7.2 Partial autocorrelation 5 year lag

for (v in unique(test$vatnsfall)){
  #acf(na.omit(test$m2[test$vatnsfall==v]), main=v)
  pacf(na.omit(test$mLm5[test$vatnsfall==v]), main=v)
}

8 Results

  1. What probability distribution should we use?
    • log-Normal distribution seems to adequate. Might be worth it to take a closer look at smaller stocks
  2. How should we estimate parameters? Is the variation different between rivers?
    • We should work on log-scale, geometric mean instead of arithmetic. The standard deviation seems to similar among river (0.5). More variation for small stocks biological or a “measurement error”?
  3. Do we need to take temporal correlation into account?
    • Strong correlation with 5 year mean and 10 year mean also significant. One year auto-correlation in some rivers (2SW rivers)
    • 10 year mean used in the previous edition of the risk assessment